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ABSTRACT 

A high-order  compact-differencing  and  Pade-type  filtering  algorithm,  coupled  with  either  the  standard  fourth-order  Runge-Kutta 
scheme  or  with  a sub-iterative  implicit  method  is  developed  and  implemented  to  simulate  aeroacoustic  phenomena  on  curvilinear 
geometries.  Several  issues  pertinent  to  the  use  of  such  schemes  are  addressed.  The  impact  of  mesh  stretching  in  the  generation  of 
high-frequency  spurious  modes  is  examined  and  the  need  for  a discriminating  higher-order  spatial  filter  procedure  is  established  and 
resolved.  The  incorporation  of  these  filtering  techniques  also  permits  a robust  treatment  of  outflow  radiation  condition  by  taking 
advantage  of  energy  transfer  to  high  frequencies  caused  by  rapid  mesh  stretching.  Computations  demonstrate  that  these  algorithmic 
components  are  also  suitable  for  interface  treatments  created  in  domain-decomposition  strategies.  For  three-dimensional 
computations,  special  metric  relations  are  employed  to  assure  the  fidelity  of  the  scheme  in  curvilinear  and  dynamic  meshes. 
Application  to  several  benchmark  2D  and  3D  problems  demonstrates  the  success  of  the  overall  computational  approach. 

INTRODUCTION 


An  important  aspect  of  both  civilian  and  military  aircraft  operation  is  the  impact  of  aerodynamically  generated  sound 
on  communities  and  structures.  Examples  of  applications  of  current  interest  include  weapon  cavity  acoustics,  jet 
screech,  sonic  boom,  cabin  noise,  and  sound  generated  by  blade/vortex  interactions.  The  relatively  new  field  of  time- 
domain  computational  aeroacoustics  (CAA)  focuses  on  the  accurate  prediction  of  aerodynamic  sound  generated  by 
airframe  components  and  propulsion  systems,  as  well  as  on  its  propagation  and  far-field  characteristics.  Both  aspects  of 
the  problem  (i.e.  sound  generation  and  propagation)  are  extremely  demanding  from  a time-domain  computational 
standpoint  due  to  the  large  number  of  grid  points  and  small  time-steps  that  are  typically  required.  Therefore,  for 
realistic  aeroacoustic  simulations  to  become  more  feasible,  higher-order  accurate  and  optimized  numerical  schemes  are 
sought  in  order  to  reduce  the  required  number  of  grid  points  per  wavelength  while  still  ensuring  tolerable  levels  of 
numerically-induced  dissipation  and  dispersion.  These  strict  simulation  requirements  are  shared  by  other  time-domain 
wave  propagation  applications  such  as  computational  electromagnetics. 

Recent  reviews  of  computational  aeroacoustics  have  been  given  by  Tam[l]  and  Wells  and  Renaut[2]  who  discuss 
various  numerical  schemes  currently  popular  in  CAA.  These  include  among  others  the  dispersion-relation-preserving 
(DRP)  scheme  of  Tam  and  Webb[3],  the  family  of  high-order  compact  differencing  schemes  of  Lele[4]  and  Essentially 
Non-Oscillatory  (ENO)  schemes[5].  The  DRP  and  compact  schemes  are  centered  non-dissipative  schemes,  a property 
which  is  desirable  for  linear  wave  propagation.  However,  their  inherent  lack  of  numerical  dissipation  results  in 
spurious  numerical  oscillations  and  instability  in  practical  applications  involving  general  geometries,  approximate 
boundary  conditions  or  non-linear  features.  In  the  DRP  approach  for  instance,  artificial  selective  damping  must  be 
employed  under  these  conditions!!].  While  quite  robust,  standard  upwind  and  upwind-biased  formulations  are 
undesirable  for  situations  involving  linear  wave  propagation  due  to  their  excessive  dissipation.  To  overcome  this 
difficulty,  higher-order  upwind  or  ENO  approaches  must  be  employed.  The  above  spatial  semi-discretizations  are 
typically  combined  with  high-order  explicit  time-integration  methods  such  as  the  multi-stage  Runge-Kutta  procedure. 
In  addition  to  the  spatial  and  temporal  discretizations,  another  critical  aspect  in  CAA  simulations  is  the  accurate 
treatment  of  the  physical  and  computational  boundary  conditions.  A recent  review  of  radiation,  outflow  and  wall 
boundary  treatments  has  been  provided  in  Ref.  6. 

This  paper  focuses  on  the  development  and  evaluation  of  a high-order  computational  methodology  for  aeroacoustic 
simulation  on  general  geometries.  There  are  two  primary  components  to  the  algorithm  chosen  in  the  present  work.  The 


Paper  presented  at  the  RTO  A VT  Symposium  on  "Ageing  Mechanisms  and  Control: 
Part  A -Developments  in  Computational  Aero-  and  Hydro-Acoustics", 
held  in  Manchester,  UK,  8-11  October  2001,  and  published  in  RTO-MP-079(1). 


{SYA)25-2 


first  is  the  spatial  differencing  scheme,  for  which  the  choice  rests  primarily  on  Padc-typc  tridiagonal-based  fourth-  and 
sixth-order  compact-difference  formuIas[4].  As  compared  to  explicit  schemes,  this  approach  incurs  an  increase  in 
computational  cost  but  lowers  the  dispersion  error  and  reduces  the  number  of  points  near  the  boundary  where  special 
formulations  are  required.  As  noted  earlier,  the  non-dissipative  nature  of  centered  schemes  makes  them  susceptible  to 
the  unrestricted  growth  of  spurious  perturbations.  This  issue  is  resolved  through  the  incorporation  of  high-order  spatial 
low-pass  filters.  The  filters  are  based  on  Pade-type  formulas  requiring  the  solution  of  tridiagonal  systems  of  equations. 
The  expressions  are  taken  from  Ref.  7,  where  up  to  10!tl-order  filters  were  employed  to  stabilize  finite-volume 
electromagnetics  calculations,  and  have  been  successfully  extended  to  solve  the  Navier-Stokes  equations[8].  Time 
integration  is  achieved  through  the  use  of  a either  the  standard  4th-order  Runge-Kutta  algorithm  or  an  implicit  sub- 
iterative  third-order  method.  The  solver  is  written  in  general  curvilinear  coordinates  in  order  to  handle  non-trivial 
geometries.  The  various  components  of  the  numerical  procedure  are  described  below  with  particular  emphasis  on 
metric  evaluation  for  3-D  and  dynamic  meshes,  multi-domain  interface  treatment  and  the  extension  of  the  filtering 
scheme  to  treat  outflow  radiation  boundaries.  The  accuracy  properties  of  the  scheme,  and  the  efficacy  of  the  high-order 
low-pass  filter  approach,  are  highlighted  by  considering  several  benchmark  test-cases  on  single  and  multiple  domains. 

METHODOLOGY 


GOVERNING  EQUATIONS 

In  order  to  develop  a procedure  suitable  for  nonlinear  fluid  dynamic,  aeroacoustic  and  aeroelastic  applications  over 
complex  configurations,  the  full  Navier-Stokes  equations  are  selected  and  are  cast  in  strong  conservative  form 
introducing  a general  time-dependent  curvilinear  coordinate  transformation  (x,y,z,t)  ->  (£,  7}  In  vector  notation, 
and  employing  non-dimensional  variables,  these  equations  can  be  written  as: 

d[u/j)/dr  + dF/dZ  + c)G/dri  + dH/dC  = (l/Re)[dFv/dZ+dGv/dTi+dHv/dC]+S/J  0) 


where  F,  G,  H,  <7„  Hv  are  the  fluxes,  5 denotes  an  imposed  acoustic  source,  and  J is  the  coordinate  transformation 
Jacobian,  which  for  dynamic  meshes,  is  also  a function  of  time.  The  specific  form  the  fluxes  appearing  in  Eq.  (1)  can 
be  found  for  instance  in  Ref.  9. 


SPATIAL  DISCRETIZATION 

A finite-difference  approach  is  employed  to  discretize  the  above  equations.  This  choice  is  motivated  by  the  relative 
ease  of  formal  extension  to  higher-order  accuracy.  For  any  scalar  quantity,  (j) , such  as  a metric,  flux  component  or 


flow  variable,  the  spatial  derivative  0 is  obtained  in  the  transformed  plane  by  solving  the  tridiagonal  system: 


+ <j>  i +r<p  n\  — b 


-+  a 


0M-&- 


(2) 


where  r,  a and  b determine  the  spatial  properties  of  the  algorithm.  The  formula  yields  the  compact  five-point,  sixth- 
order  (C6)  and  three-point  fourth-order  (C4)  schemes  with  r= 1/3,  o=14/9,  £=1/9  and  7=1/4,  n=3/2,  £= 0 respectively. 
Equation  (2)  also  incorporates  the  standard  explicit  fourth-order  (E4)  and  second-order  (E2)  approaches  for  which  the 
coefficients  are  (T=  0,  cr=4/3,  £=-1/3)  and  (/=  0,  a=l,  £= 0)  respectively.  The  dispersion  characteristics  and  truncation 
error  of  the  above  schemes  can  be  found  in  Refs.  4 and  10.  It  should  be  noted  that  for  a given  order  of  accuracy,  the 
compact  schemes  are  significantly  superior  to  their  explicit  (non-compact)  counterparts. 


FILTERING  SCHEME 

Compact-difference  discretizations,  like  other  centered  schemes,  are  non-dissipative,  and  therefore  susceptible  to 
numerical  instabilities  due  to  the  unrestricted  growth  pf  high-frequency  modes.  These  difficulties  originate  from  several 
sources  including  mesh  non-uniformity,  approximate  boundary  conditions  and  nonlinear  flow  features.  In  order  to 
extend  the  present  solver  to  practical  simulations,  while  retaining  the  improved  accuracy  of  the  spatial  compact 
discretization,  a high-order  implicit  filtering  technique  [7,8]  is  incorporated.  If  a typical  component  of  the  solution 

vector  is  denoted  by  <p  , filtered  values  (j)  satisfy, 

+t+a  A i = S + ) (3) 

n=0  *■ 

Equation  (3)  is  based  on  templates  proposed  in  Ref.  4 and  with  proper  choice  of  coefficients,  provides  a 2/Vth-order 
formula  on  a 2N+1  point  stencil.  The  N+]  coefficients,  au,  at,  %,  are  derived  in  terms  of  Of  with  Taylor-  and 
Fourier-series  analyses  and  are  found  in  Refs.  8,  10  and  11,  along  with  the  corresponding  spectral  filter  responses.  The 
adjustable  parameter  ot/  satisfies  the  inequality  -0.5  < a f < 0.5,  with  higher  values  of  ct/  corresponding  to  a less 
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dissipative  filter.  In  multi-dimensional  problems,  the  filter  operator  is  applied  sequentially  in  each  coordinate  direction. 
The  specified  interior  filter  formula  is  denoted  by  appending  its  designation  to  that  of  the  scheme.  For  example,  C6F10 
designates  the  sixth-order  compact  scheme  combined  with  a tenth-order  filter.  For  the  near-boundary  points,  the 
filtering  strategies  described  in  Refs.  8 and  1 1 are  employed. 

METRIC  EVALUATION 

The  extension  of  high-order  schemes  to  non-trivial  3-D  geometries  demands  that  issues  of  ffeestream  preservation  and 
metric  cancellation  be  carefully  addressed.  These  errors,  which  arise  in  finite-difference  discretizations  of  governing 
equations  written  in  strong-conservation  form,  can  catastrophically  degrade  the  fidelity  of  standard  second-order  as 
well  as  higher-order  approaches[8,12].  In  deriving  the  flow  equations  in  strong-conservation  form,  the  following 
metric  identities  have  been  implicitly  invoked, 

h = (fi  /j)(  + {l)xlJ)rj  + (C x/J)i ; = 0 

h = tey/  -Of  + + (Cy/ J\  = 0 (4) 

h = (UJk  + MJ)n  + (UJ\  = o 
h - (1/ J)r  + (tt/Jk  + + ( Ci/Jk  = 0 

In  Ref.  8 it  was  shown  that  on  highly  distorted  (static)  curvilinear  2-D  meshes,  the  compact  scheme  exhibits  ffeestream 
preservation  when  the  metrics  are  evaluated  with  the  same  finite-difference  expressions  as  those  employed  for  the 
fluxes.  It  was  also  demonstrated  that  the  practice  of  prescribing  analytic  metrics  on  stretched  curvilinear  grids  can  lead 
to  unacceptable  errors  and  therefore  should  in  general  be  avoided.  The  previous  straightforward  approach  of  calculating 
the  metrics,  although  effective  in  2-D,  fails  to  provide  metric  cancellation  for  general  3-D  curvilinear  configurations. 
To  illustrate  this  point,  consider  the  standard  metric  relations: 

Cx/J .=  - y^Zff 

Vx/J  = (5) 

Cx/J  = V£zti  — yr)z£ 

associated  with  the  identity  I,.  Evaluation  of  the  y and  z derivatives  in  the  previous  expressions  using  explicit  or 
compact  centered  schemes  does  not  satisfy  I,,  and  as  a result,  significant  grid-induced  errors  may  appear[12,13].  In 
order  to  extend  the  high-order  compact  scheme  to  general  geometries,  the  metric  terms  are  rewritten  prior  to 
discretization  in  the  equivalent  (conservative)  form[14]: 

Cx/J  = {vnz)(  - {y^h 

Vx/J  = ~ (6) 

Cx/J  = {y^)v  ~ (yyzk 

Similar  expressions  are  employed  for  the  remaining  metric  terms  in  order  to  satisfy  the  identities  I2  and  I3  above.  When 
the  transformation  metrics  are  recast  in  this  manner,  and  the  derivatives  are  evaluated  with  the  same  high-order 
formulas  employed  for  the  fluxes,  freestream  preservation  is  again  recovered  in  general  time-invariant  3-D  curvilinear 
geometries[12]. 

For  deforming  and  moving  meshes,  the  identity  I4  must  be  also  satisfied  to  eliminate  metric  cancellation  errors  and  to 
ensure  ffeestream  preservation.  This  metric  identity  is  referred  to  in  the  literature[14]"as  the  Geometric  Conservation 
Law  ( GCL ).  For  the  time-integration  methods  employed  in  this  work,  the  time-derivative  term  in  Eq.  (1)  is  split  using 
chain-rule  differentiation  as  follows: 

(U/J)T  = {1/J)UT  + U(l/J)T  (7) 

Rather  than  attempting  to  compute  the  time  derivative  of  the  inverse  Jacobian  directly  from  the  grid  coordinates  at 
various  time  levels  (either  analytically  or  numerically),  we  simply  invoke  the  GCL  identity  I4  to  evaluate  (1/J)r,  i.e. 

(i/J)r  = -m/J^  + {rk/J)n  + (Ct/J\\ 

where 

Ct/J  = -M£r  /J)  + yACy/J)  + ZriCz/J)} 

m/J  = -W<Vx/J)  + VriVy/J)  + zT{Vz/J)\ 

Ct/J  = -WkCxU)  + yr(Cy/J)  + Zr(Cz/J)\ 

For  the  case  of  an  analytically  prescribed  dynamic  mesh  transformation,  the  grid  speeds  (xT,yt,zT)  are  obtained  from  the 
corresponding  analytic  expressions.  An  example  in  which  the  grid  speeds  are  known  analytically  corresponds  to  the 


(8) 

(9) 
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case  of  a maneuvering  wing  when  the  entire  mesh  is  rotated  in  a rigid  fashion.  In  many  practical  applications  involving 
deforming  meshes  (e.g.  dynamic  aeroelastic  simulations),  the  grid  speeds  are  not  known  a priori,  and  must  therefore  be 
approximated  to  the  desired  degree  of  accuracy  employing  the  evolving  grid  coordinates  at  several  time  levels.  As 
demonstrated  in  Ref.  15,  the  high-order  method  retains  its  superior  accuracy  on  rapidly  distorting  meshes  when  the 
procedure  outlined  above  is  incorporated  for  the  time  metrics. 


TIME-INTEGRATION  SCHEME 

Two  different  time-integration  approaches  are  incorporated  in  the  present  family  of  solvers.  For  wave  propagation 
applications,  the  equations  are  integrated  in  time  with  the  classical  fourth-order  four-stage  Runge-Kutta  method  (RK4). 
The  scheme  is  implemented  in  low  storage  form  requiring  three  levels  of  storage.  For  highly  stretched  grids  exhibiting 
very  small  volumes,  the  stability  constraint  of  the  explicit  time-marching  scheme  is  found  to  render  the  approach  too 
restrictive  and  inefficient.  Therefore,  the  implicit,  approximately-factored  method  of  Beam  and  Warming[16]  is  also 
incorporated  and  augmented  through  the  use  of  Newton-like  subiterations  in  order  to  achieve  high-order  time  accuracy. 
In  delta  form,  the  second-order  form  of  the  scheme  (denoted  as  BW2)  may  be  written  as: 


b (x)< * (sj- - m )]• (f  )■ •>{%■ - m )]> « b (x)' * ($- - m )]- * 

- - (x)  [(sa)  (W  • - ■ + « «--)  +A  (p  - in) + «,  (<? - id)  + (p  - in) + *] 


(10) 


The  spatial  derivatives  in  the  implicit  operators  are  represented  using  standard  second-order  centered  approximations 
whereas  high-order  discretizations  are  employed  for  the  explicit  side.  Nonlinear  artificial  dissipation  terms[17]  are  also 
appended  to  the  implicit  operator  to  enhance  stability.  In  addition,  for  improved  efficiency,  the  approximately-factored 
scheme  is  recast  in  diagonalized  form[18].  Any  degradation  in  solution  accuracy  caused  by  the  spatial  second-order 
implicit  operators,  artificial  dissipation  and  the  diagonal  form  is  eliminated  through  the  use  of  sub-iterations  (typically, 
three  sub-iterations  are  used  per  time  step).  By  changing  the  number  of  time  levels  employed  to  evaluate  the  time- 
derivative  term  appearing  in  the  RHS  of  Eq.  (10),  first  (BW1)  and  third  {BW3)  order  accurate  forms  of  the  implicit 
algorithm  can  be  constructed. 


MULTI-DOMAIN  STRATEGY 

Domain-decomposition  techniques  constitute  an  important  component  of  modem  computational  strategy.  Due  to  their 
spatially  implicit  nature,  Pade-type  schemes  are  more  difficult  to  utilize  in  a multi-domain  environment  than  explicit 
methods.  However,  a finite-size  overlap  can  be  employed  with  the  present  compact/filtering  methodology  to  generate  a 
powerful  approach  applicable  to  complicated  curvilinear  meshesfl  1,12],  Figure  1 depicts  schematically  the  case  of  a 
simulation  employing  two  sub-domains.  The  original  domain  of  computation  is  divided  into  two  parts  to  be  distributed 
to  different  processors.  Each  sub-domain  is  supplemented  with  several  points  from  the  adjacent  sub-domain  to  form  an 
overlap  region,  whose  details  for  a five-point  vertical  overlap  are  also  shown  in  Fig.  1 . Although  the  overlap  points  are 
collocated,  they  have  been  shown  slightly  staggered  for  clarity.  Each  vertical  line  is  denoted  by  its  /-index.  Data  is 
exchanged  between  adjacent  sub-domains  at  the  end  of  each  sub-iteration  of  the  implicit  scheme  (or  each  stage  of 
RK4 ),  as  well  as  after  each  application  of  the  filter.  The  values  at  points  1 and  2 of  Mesh  2 are  set  to  be  identically 
equal  to  the  corresponding  updated  values  at  points  IL-4  and  IL-3  of  Mesh  1 . Similarly,  reciprocal  information  is 
transferred  through  points  4 and  5 of  Mesh  2 which  “donate”  values  to  points  IL-1  and  IL  of  Mesh  1.  More  details  on 
the  accuracy  and  robustness  of  the  presentmulti-domain  approach  can  be  found  in  Ref.  11. 


RESULTS 


IMPACT  OF  SPATIAL  FILTERING  ON  STRETCHED  MESHES 

In  practical  scattering  simulations,  grid  stretching  is  usually  employed  in  order  to  reduce  required  computational 
resources,  as  well  as  to  permit  the  use  of  approximate  farfield  boundary  conditions.  Therefore,  the  performance  of  high- 
order  schemes  on  general  stretched  meshes  must  be  carefully  examined. 

An  enlightening  analysis  of  the  behavior  of  a smooth  solution  as  it  passes  through  a sudden  mesh  coarsening  has  been 
presented  by  Vichnevetsky[19]  for  the  1-D  advection  equation  semi-discretized  with  the  standard  second-order 
centered  scheme.  This  analysis  indicates  that  although  total  energy  is  preserved,  at  the  grid  coarsening  interface,  a 
significant  portion  of  the  energy  is  deposited  on  a reflected  solution  composed  primarily  of  odd-even  modes  and 
modulated  by  a smooth  envelope.  This  reflected  energy  propagates  upstream  (i.e.  with  negative  group  velocity)  and  in 
most  circumstances  if  left  unchecked  has  the  potential  of  ultimately  contaminating  the  genuine  solution.  However, 
since  reflections  due  to  grid  stretching  are  characterized  by  high-frequency  modes,  they  can  be  easily  removed  by  the 
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high-order  low-pass  sparial  filter.  The  effectiveness  of  the  low-pass  filter  in  controlling  spurious  reflections  without 
degrading  the  fidelity  of  the  solution  was  demonstrated  previously  in  Ref.  12. 

The  combination  of  grid  stretching  with  a discriminating  low-pass  filter  may  be  exploited  as  an  alternative  procedure 
for  outflow  boundary  treatment.  By  employing  a large  rate  of  stretching,  a significant  amount  of  energy  can  be 
reflected  at  the  grid  coarsening  interfaces.  Provided  this  reflected  energy  is  deposited  into  high-frequency  modes  (in  the 
fine  mesh  region),  it  can  then  be  subsequently  eliminated  by  the  baseline  low-pass  filter  without  contaminating  the 
genuine  solution.  Furthermore,  the  transmitted  energy  is  also  quickly  dissipated  by  the  high-order  filter,  as  the 
transmitted  solution  features  are  represented  by  a diminishing  number  of  points  per  wave  in  the  stretching  mesh.  This 
proposed  method  eliminates  the  need  for  more  sophisticated  (and  perhaps  more  restrictive)  boundary  conditions  at  the 
expense  of  extending  the  computational  domain.  Although  the  proposed  approach  has  a mathematical  foundation  (at 
least  based  on  the  1-D  analysis  of  Ref.  19),  its  implementation  is  highly  empirical,  and  therefore  its  utility  must  be 
evaluated  in  the  context  of  practical  applications. 

As  a severe  test  case,  consider  the  propagation  of  acoustic  waves  in  the  grid  of  Fig.  2a.  This  mesh  is  uniform 
(Ax=Ay=0.05)  in  the  center  of  the  computational  domain  (- 3 < x,y  < 3).  Outside  of  this  resolved  region,  Ax  and  Ay  are 
increased  abruptly.  An  acoustic  source  is  specified  according  to  the  expression 

S(x,y,t)  = exp{(ln2)  [(x-xj2+(y-yj2/ b2]}  sin(a>t)  f(t)  (1 1) 

f(t)  =min{  1.0,  (t/t0)3} 

with  xc=  yc=  0,  tu  = 5k,  b = 0.2  and  t„  = 4.  A snapshot  of  the  pressure  is  displayed  in  Fig.  2b.  It  is  apparent  that  the 
acoustic  energy  reflected  at  the  grid-coarsening  interface  is  almost  completely  annihilated  by  the  high-order  filter.  The 
corresponding  plot  of  the  instantaneous  pressure  along  the  diagonal  (x  - y)  is  shown  in  Fig.  2c  and  indicates  that  the 
transmitted  energy'  diffuses  rapidly  in  the  coarse  mesh  region.  For  this  square  grid,  the  circular  waves  cross  the  grid 
interface  at  a varying  angle  of  incidence  without  apparent  anisotropies  being  introduced.  More  quantitative  tests 
demonstrating  the  suitability  of  the  present  far  field  radiation  treatment  can  be  found  in  Ref.  13  for  a transient  acoustic 
pulse  and  a convective  vortical  disturbance. 

SCATTERING  OF  AN  ACOUSTIC  PULSE 

In  order  to  validate  the  present  approach  for  curvilinear  geometries,  we  select  as  a test  case  the  benchmark  problem 
denoted  as  Category  I,  Problem  2 in  the  2nd  CAA  Workshop  of  Ref.  20.  This  configuration  (Fig.  3)  describes  the 
scattering  off  a circular  cylinder  of  a prescribed  initial  pressure  pulse.  The  pulse  at  t = 0 is  given  by  the  expression 

p =p0  + eexp{(ln2)  [(x-xf+(y-yj2/b2]}  (12) 

with  xc=4,  yc  =0,  e=0.01,  b=0.2. 

Along  the  cylinder  surface,  standard  inviscid  wall  boundary  conditions  are  employed.  Since  the  configuration  is 
symmetric,  only  the  upper  half  of  the  domain  is  considered,  and  symmetry  conditions  are  invoked  along  6 = (f.  J80°. 
As  indicated  in  Fig.  3a,  the  grid  is  stretched  for  r/D  > 7 using  a constant  stretching  factor  of  1.1.  Since  in  this  coarse- 
mesh  region  the  outgoing  pulse  isclissipated  by  the  baseline  filter,  a simple  extrapolation  condition  is  suitable  along  the  - 
far  field  boundary.  A polar-type  grid  of  size  156X175  is  employed,  and  the  solution  is  advanced  in  time  with  a non- 
dimensional  time  step  At  = 0.004. 

Pressure  contours  depicting  pulse  propagation  and  reflection  off  the  cylinder  surface  are  shown  in  Fig.  3 at  several 
instants  in  time  for  the  C6F10-RK4  scheme.  The  history  of  the  pressure  at  a selected  point  (r/D= 5,  0=9(f  ) and  for 
several  computations  is  presented  in  Fig.  4.  On  this  relatively  coarse  mesh,  the  C6F10-RK4  and  C4F8-RK4  schemes 
with  a filter  coefficient  oy=  0.4  are  observed  to  be  in  good  agreement  with  the  exact  solution[20].  However,  on  this 
mesh,  the  C4  scheme  combined  with  an  8lh-order  explicit  filter  ( i.e . af  = 0 ) displays  appreciable  error.  This  highlights 
the  improved  accuracy  of  the  implicit  filter  formulation. 

Results  computed  using  the  implicit  sub-iterative  time-marching  scheme  are  provided  in  Fig.  5.  These  results  were 
obtained  with  a non-dimensional  time  step  At  = 0.004  and  employing  two  sub-iterations.  The  first  order  time 
advancement  method  (BW1)  is  found  to  be  highly  dissipative.  However,  significant  improvement  is  achieved  by 
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employing  the  second-order  implicit  approach  ( BW2 ).  Finally,  the  third-order  iterative  implicit  method  (B W3)  is  in  very 
close  agreement  with  the  fourth-order  explicit  (RK4)  solution  (see  figure  inset). 

In  order  to  demonstrate  the  advantages  of  the  implicit  time-marching  scheme,  we  consider  the  previous  pulse  scattering 
example  for  the  case  of  an  elliptic  cylinder  of  aspect  ratio  10  (Fig.  6).  As  Fig.  6a  indicates,  the  grid  is  highly  clustered 
near  the  small  radius  of  curvature  at  the  edge  of  the  ellipse.  This  results  in  small  grid  cells  and  consequently  in  a 
significant  reduction  of  the  allowable  time  step  for  the  explicit  method.  Calculations  employing  the  RK4  scheme  were 
found  to  be  unstable  even  for  a non-dimensional  time  step  as  low  as  At  = 0.000125  which  correspond  to  a maximum 
CFL  number  of  approximately  0.96.  Therefore  it  becomes  apparent  that  the  explicit  time-marching  method  is  very 
inefficient  for  general  geometries  wherein  highly  clustered  meshes  are  necessary  (e.g.  sharp  edges).  By  contrast,  the 
third-order  implicit  iterative  scheme  was  found  to  remain  stable  for  the  baseline  time  step  (At  = 0.004)  which 
corresponds  to  CFLmai  ~ 31.  The  corresponding  instantaneous  pressure  contours  obtained  with  the  implicit  method  at  t 
= 1.5  are  displayed  in  Fig.  6b. 

MULTI-DOMAIN  SCATTERING  SIMULATION 

This  case  corresponds  to  Category  I,  Problem  1 in  Ref.  20,  and  considers  the  scattering  from  a circular  cylinder  of  a 
periodic  acoustic  source  (Fig.  7a).  This  example  constitutes  a more  stringent  test  of  the  algorithm  and  boundaiy 
conditions  since  an  asymptotic  periodic  solution  must  be  attained,  and  long-term  numerical  stability  demonstrated.  In 
addition,  in  order  to  highlight  the  capability  of  the  present  numerical  approach  to  treat  a multiple-domain  situation,  the 
single-domain  grid,  consisting  of  361  X 321  points,  is  split  along  6 - 90°,  where  extra  £-lines  are  added  to  form  a five- 
point  overlap  as  in  Fig.  1 . Solutions  are  advanced  separately  on  each  sub-domain,  and  information  is  exchanged  at  the 
overlap  points  in  the  manner  previously  discussed.  The  C6  scheme  is  employed  for  interior  points  along  with  fourth- 
and  fifth-order  compact  operators  at  the  boundary  and  next-to-boundary  points  respectively  whereas  RK4  is  utilized  for 
time-integration.  In  the  interior  of  each  sub-domain,  a 1 0lh-order  filter  is  utilized  whereas  high-order  one-sided 
techniques,  are  invoked  near  boundaries.  For  all  filter  operators,  the  coefficient  aj~  0.45  is  specified. 

Figure  7b  displays  instantaneous  pressure  contours  in  the  vicinity  of  the  cylinder.  It  is  apparent  that  the  pressure  waves 
cross  the  grid  interface  without  producing  any  noticeable  disruptions  of  the  interference  pattern  even  though  pressure 
waves  generated  by  the  source  propagate  through  the  overlap  region  in  an  oblique  direction  to  the  zonal  interface.  A 
quantitative  comparison  of  the  single-domain,  multiple-domain  and  analytic  solutions  is  given  in  Fig.  7c  in  terms  of  the 
directivity  of  the  radiated  sound  at  r/D  = 5.  The  directivity  obtained  with  a 6lh-order  near-boundary  filter  in  the  overlap 
region  is  essentially  the  same  as  the  corresponding  single-domain  baseline  solution,  and  both  results  are  in  excellent 
agreement  with  the  theoretical  solution.  These  results  highlight  the  potential  of  the  present  high-order  methodology  for 
domain-decomposition  applications  on  parallel  computers. 

ACOUSTIC  PROPAGATION  ON  3-D  CURVILINEAR  MESHES 

This  validation  case  considers  the  propagation  of  a 3-D  spherical  pulse  in  a curvilinear  mesh  in  order  to  examine  metric 
cancellation  errors.  A three-dimensional  curvilinear  61X61X61  mesh  is  generated  using  the  equations  given  in  Ref.  13, 
providing  a grid  in  which  the  metric  identities  (Eq.  (4))  are  not  trivially  satisfied.  A constant  i^-plane  of  the  mesh  at  the 
location  of  maximum  deformation  is  shown  in  Fig.  8a.  The  imposed  grid  undulations  are  resolved  with  approximately 
1 5 points  per  wave.  An  initial  pressure  pulse  is  prescribed  by 

p~po  + eexpf  (-ln2)  (x ?+y2+z2)/9) 

where  e-  0.01. 

The  pulse  propagation  problem  is  computed  with  a At  = 0.004  using  the  C4F10-RK4  algorithm  with  (Xf  = 0.49.  Figure 
8b  displays  the  calculated  pressure  contours  on  a plane  through  the  center  of  the  spherical  pulse  at  t=10  for  the  case 
when  the  metrics  are  evaluated  in  the  standard  fashion  of  Eq.  (5).  It  is  apparent  that  significant  distortion  of  the  wave 
front  occurs  due  to  the  lack  of  freestream  preservation  ( i.e . metric  cancellation  errors).  The  perturbation  pressure  along 
the  grid  line  i-j-31,  (Fig.  8d)  indicates  gross  departure  from  the  theoretical  solution.  The  results  obtained  with  the 
metric  evaluation  procedure  of  Eq.  (6)  exhibit  no  distortions  of  the  spherical  front  (Fig.  8c)  and  are  in  excellent 
agreement  with  the  exact  solution  (Fig.  8d).  Calculations  with  C6F10  (not  shown)  displayed  reduced  sensitivity  to  the 
choice  of  metric  evaluation  procedure.  This  is  in  agreement  with  the  results  of  Ref.  12,  wherein  metric  cancellation 
errors  were  shown  to  decrease  with  increasing  order  of  accuracy.  Nonetheless,  without  the  incorporation  of  Eq.  (6)  in 
the  calculation  of  the  metric  relations,  all  solutions  on  this  distorted  mesh  were  found  to  be  of  poor  quality.  Hence,  even 
for  relatively  benign  3-D  curvilinear  grids,  freestream  preservation  errors  may  swamp  acoustic  pressure  levels  unless 
proper  metric  evaluation  procedures  are  employed.  It  should  be  noted  that  analytic  metric  evaluation  (even  if  available) 
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does  not  remedy  the  situation  [8,12].  The  present  results  clearly  demonstrate  that  the  superior  performance  of  the  high- 
order  method  can  be  extended  to  the  realm  of  general  curvilinear  grids  including  moving/deforming  meshes  as 
highlighted  in  the  fluid/structure  interaction  described  next. 

ACOUSTO-STRUCTURAL-FLUID  INTERACTION  OVER  A FLEXIBLE  PANEL 

As  a final  example  of  a simulation  of  multi-disciplinary  physics  with  the  present  methodology,  consider  a transitional 
boundary-layer  flow  over  a flexible  finite  panel  embedded  in  a rigid  surface  as  shown  schematically  in  Fig.  9a.  This 
problem  is  closely  related  to  classic  panel  flutter  phenomena,  as  well  as  to  viscous  flow  over  compliant  surfaces.  The 
panel  of  length  a and  thickness  h extends  over  the  region  0.5  < x/a  < 1.5.  The  leading-edge  region  of  the  plate  (-0.5  < 
x/a  < 0.0 ) is  formed  by  an  ellipse  of  half-thickness  0.05h  ( i.e . aspect  ratio  10).  An  additional  challenge  posed  by  this 
aeroelastic  simulation  is  the  need  to  accomodate  the  surface  deflection  with  a dynamically  deforming  mesh.  The 
problem  has  been  examined  in  great  detail  in  Ref.  21,  which  should  be  consulted  for  details  regarding  boundary 
condition  implementation  and  mesh  resolution  studies.  For  brevity,  only  select  results  obtained  at  a ffeestream  Mach 
number  M = 0.8  and  Rea  = 10s  are  summarized  here  to  highlight  the  ability  of  the  method  to  capture  the  complicated 
unsteady  phenomena  under  the  influence  of  flow-induced  surface  deformation.  This  case  was  computed  employing  the 
C6F10-BW2  scheme. 

At  low  values  of  the  dynamic  pressure,  a steady  flow  is  obtained  despite  the  adverse  pressure  gradient  induced  by  the 
downward  deflection  of  the  panel.  At  higher  dynamic  pressures,  however,  a traveling-wave-flutter  phenomenon  is 
observed  as  summarized  in  Fig.  9b,  c.  The  instantaneous  panel  shapes  (not  shown)  display  a seventh-mode  oscillation 
with  a dominant  non-dimensional  frequency  St  - fa/U  = 1.52  which  is  substantially  higher  than  the  fundamental 
frequency  of  the  elastic  plate.  The  high-mode  flexural  deflections  are  observed  to  travel  along  the  panel  and  to  reflect  at 
the  panel  edges.  These  high-frequency  fluctuations  result  in  a pronounced  acoustic  radiation  pattern  above  the  vibrating 
plate,  shown  in  Fig.  9b  in  terms  of  the  instantaneous  pressure  field.  Downstream  of  the  flexible  surface,  a regular  train 
of  vortical  disturbances  is  observed  (Fig.  9c)  with  characteristic  wavelength  and  frequency  compatible  with  those  of 
Tollmien-Schlichting  (T-S)  instability.  The  traveling  wave  flutter  appeared  to  originate  from  the  coupling  of  the  T-S 
waves  with  the  panel  high-mode  transverse  fluctuations,  and  this  convective  instability  ceases  below  a critical  value  of 
Reynolds  number. 

SUMMARY 

A high-order  finite-difference  method  has  been  adapted  to  simulate  aeroacoustic  phenomena  on  curvilinear  geometries. 
The  scheme  is  based  on  4th-  and  6,h-order  compact  differencing  formulas  coupled  with  up  to  10lh-order  Pade-type  low- 
pass  spatial  filters.  Both  explicit  and  iterative  implicit  time-marching  schemes  are  considered.  The  high-order  spatial 
filtering  procedure  is  required  in  order  to  maintain  stability  in  the  presence  of  mesh  stretching  and  approximate 
boundary  conditions  while  ensuring  high  fidelity  for  wave  propagation.  The  incorporation  of  these  filtering  techniques 
also  permit  a robust  treatment  of  outflow  radiation  condition  by  taking  advantage  of  energy  transfer  to  high-frequencies 
caused  by  rapid  mesh  stretching.  The  potential  of  the  procedure  for  parallel  implementation  is  demonstrated  by 
successful  application  to  multi-domain  scattering  cases.  Special  metric  procedures  are  also  shown  to  be  essential  in  the 
simulation  of  acoustic  phenomena  in  general  curvilinear  three-dimensional  meshes. 
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Fig.  1.  Schematic  of  overlap  strategy  in  high- 
order  domain  decomposition. 


Fig.  2.  Computation  of  2-D  acoustic  source  on  mesh  with  abrupt  stretching,  (a)  Grid, 
(b)  Pressure  contours,  (c)  Pressure  along  diagonal. 
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Fig.  3.  Pressure  contours  at  various  instants  for  acoustic 
pulse  scattered  by  a circular  cylinder. 


Fig.  4.  Computed  pressure  history  at  select  Fig.  5.  Comparison  of  implicit  and 

point  for  several  spatial  discretizations  and  explicit  time-marching  schemes, 

the  RK4  scheme. 
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Fig.  6.  Acoustic  pulse  scattered  by  elliptic  cylinder,  (a)  Mesh  , (b)  Instantaneous  pressure 
contours  at  t = 1.5  fox~C6F10-BW3  scheme. 


Fig.  7.  Multi-domain  scattering  simulation,  (a)  Configuration,  (b)  pressure  contours  near  mesfi 
overlap,  (c)  Directivity  at  r/D  - 5. 
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Fig.  8.  Propagation  of  spherical  pulse  on  a curvilinear  mesh, 
(a)  Grid,  (b)  and  (c)  Pressure  contours  for  standard  and  new 
metrics  respectively,  (d)  pressure  distribution  along  a grid 
line  through  the  center  of  the  pulse  at  t -10. 
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Fig.  9.  Transitional  boundary-layer  flow  over  a fluttering  elastic  panel,  (a)  Configuration,  (b) 
Instantaneous  pressure,  (c)  Vorticity. 


